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& ; Abstract 

We present a general computational scheme based on molecular dynamics (m.d.) simulation for calculating the chemical 
potential of adsorbed molecules in thermal equilibrium on the surface of a material. The scheme is based on the calculation of 
the mean force in m.d. simulations in which the height of a chosen molecule above the surface is constrained, and subsequent 
""j" 1 | integration of the mean force to obtain the potential of mean force and hence the chemical potential. The scheme is valid at 

■ any coverage and temperature, so that in principle it allows the calculation of the chemical potential as a function of coverage 
and temperature. It avoids all statistical mechanical approximations, except for the use of classical statistical mechanics for 

O , the nuclei, and assumes nothing in advance about the adsorption sites. From the chemical potential, the absolute desorption 

O ■ rate of the molecules can be computed, provided the equilibration rate on the surface is faster than the desorption rate. We 
apply the theory by ab initio m.d. simulation to the case of H2O on MgO (001) in the low-coverage limit, using the Perdew- 

I , Burke-Ernzerhof (PBE) form of exchange-correlation. The calculations yield an ab initio value of the Polanyi-Wigner frequency 

■ prefactor, which is more than two orders of magnitude greater than the value of 10 13 s _1 often assumed in the past. Provisional 
' comparison with experiment suggests that the PBE adsorption energy may be too low, but the extension of the calculations to 
, higher coverages is needed before firm conclusions can be drawn. The possibility of including quantum nuclear effects by using 
' path-integral simulations is noted. 

q ■ 1 Introduction 

i> : 

, Ab initio modelling based on density functional theory (DFT) has had a major impact on the understanding of molecular 
• • . adsorption and desorption at surfaces [H [2] - However, most of this modelling has been of the static kind, in which structural 
^ ' relaxation is used to calculate the energies of chosen adsorbate geometries. For the interpretation of some types of experiment, 
, static calculations may suffice, but in many practical situations entropy effects cannot be ignored. This is clearly true when 
one considers surface phase diagrams or adsorption isotherms, but it is also true for rate processes such as thermal desorption. 
In a recent short publication [3], we showed the possibility of using DFT molecular dynamics simulation to calculate absolute 
desorption rates, with full inclusion of all entropy effects. Our aims here are to describe the general theory underlying that 
work, to outline its relevance to the calculation of adsorption isotherms, to present the simulations themselves in more detail, 
and to indicate several developments that we plan to explore in future papers. 

Although most DFT surface modelling has been static, the general idea of ab initio statistical mechanics (AISM) applied 
to condensed matter goes back well over 10 years [HE]. Over that period, methods for calculating the ab initio free energy 
of liquids and anharmonic solids have become well established 6, 7, 8 . In surface science, one of the earliest publications 
describing a form of AISM was that of Stampfl et al. 9. on the O/Ru (0001) system, in which they successfully calculated 
thermal desorption spectra, heat of adsorption and the surface phase diagram. 

All previous AISM work that we know of on surface adsorbates has been based on lattice-gas schemes, which we specifically 
wish to avoid here. To make this point clear, we recall that there are two fundamentally different approaches to the statistical 
mechanics of adsorbates. In lattice gas theories, it is assumed from the outset that the adsorbate atoms or molecules (for 
brevity, we shall simply say 'molecules') can occupy only specified surface sites. Calculations are then formulated in terms 
of the occupancies of these sites and transition rates between them, the model parameters being fixed either by appeal to 
experiment or (more recently) to DFT calculations. A completely different approach is to regard the adsorbate molecules as 
forming a two-dimensional fluid interacting with the substrate. The methods we shall describe resemble the second approach, 
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and are also close to the AISM methods developed for bulk liquids. A crucial point is that these bulk AISM methods are 
designed to deliver free energies or chemical potentials that are free of statistical-mechanical errors (more precisely, these errors 
can be systematically reduced to any desired tolerance) [101 111] , Here, we have the same aim for the AISM of adsorbates: 
the statistical mechanical errors should be systematically controllable, so that the only remaining error is due to the DFT 
exchange-correlation approximation. In a lattice-gas approach, by contrast, the adoption of a discrete set of adsorption sites 
constitutes an approximation that is difficult to eliminate, and this is why we avoid this approach. 

Coming to specifics, the modelling techniques that we wish to develop should allow the calculation of the chemical potential 
of a system of adsorbate molecules on the surface of a material at any coverage and temperature, starting from a given DFT 
exchange-correlation functional. We focus particularly on the chemical potential, since this is the key to calculating all other 
thermodynamic properties of the adsorbate, as well as the thermal desorption rate measured in temperature programmed 
desorption (TPD) experiments [12]. The techniques should not assume in advance anything about the sites occupied by the 
adsorbate molecules, or about their orientations, and should not rely on the harmonic approximation for any of the degrees of 
freedom. Ultimately, though not in this paper, we should like to be able to treat molecules that dissociate and undergo chemical 
reactions. We want to be able to include quantum nuclear effects, if necessary. Although we assume for the moment that the 
ab initio total energy function is provided by DFT, we recognise that DFT is often inaccurate, and we want the techniques to 
be generalisable to more accurate ab initio methods, such as quantum Monte Carlo [13] or high-level quantum chemistry. 

We present in the next Section a theoretical scheme which we believe will be capable of allowing many of these aims to be 
accomplished. Our strategy is to consider the system consisting of the adsorbate molecules on the surface in complete thermal 
equilibrium with the gas phase. All molecules in the system, including those in the substrate, are treated on an equal footing. 
In this situation, there is a distribution p(z), giving the probability of finding adsorbate molecules a distance z from the surface. 
This distribution can be expressed in terms of a potential of mean force (PMF): p(z) = A exp(— cj>(z)/kBT), and the PMF <f>(z) 
is the integral with respect to z of the mean force (T z ) z acting on the adsorbate molecules when they are constrained to be at 
distance z from the surface. We shall show how these simple relationships allow the calculation of the chemical potential of the 
adsorbate. We note that the same relationships also yield the absolute rate of thermal desorption as a function of coverage and 
T, provided all equilibration rates on the surface are faster than the desorption rate. We shall show that, under appropriate 
conditions, the theory yields the widely used Polanyi-Wigner equation Q3], according to which the rate of change of adsorbate 
coverage 6 (molecules per adsorption site) in a TPD experiment is given by: 

dB/dt = -f0 n exp(- AE/k B T) , (1) 

where AE and n are the activation energy and reaction order for desorption, and / is a frequency prefactor. In the analysis 
of experimental TPD data [151 116] , / is commonly treated as an empirical adjustable parameter, but the theory we present is 
capable of yielding ab initio values for /. 

The practical application of the scheme reported in this paper (Sec. [3} consists of DFT simulations of the H2O molecule 
on the perfect MgO (001) surface at several temperatures in the limit of zero coverage. We confine ourselves here to zero 
coverage, even though the theory is equally valid at all coverages, simply because the practical sampling requirements are easier 
to satisfy in this case. The target we set ourselves is to calculate the chemical potential of the adsorbate molecule in thermal 
equilibrium, with complete inclusion of entropy effects, to an accuracy of better than 50 meV for the chosen exchange-correlation 
functional. We shall see that even at low coverage entropy effects are very important, and increase the frequency prefactor for 
thermal desorption by over two orders of magnitude in the temperature range of practical interest. In Sec. [4] we discuss the 
practical problems involved in doing calculations at higher coverages, and possible ways of overcoming them; we also discuss 
the generalisation to include quantum nuclear effects, using DFT path-integral techniques. Our conclusions are summarised in 
Sec. H] 

2 Theory and techniques 

We begin our outline of theory and techniques by describing the statistical mechanical relationships that allow us to calculate 
the chemical potential of adsorbed molecules using PMF techniques. We then summarise the arguments that allow the thermal 
desorption rate to be deduced from this chemical potential in the case of fast equilibration. At the end of this Section, we 
discuss the practical issues in calculating the PMF, and we provide technical details of the DFT molecular dynamics techniques 
used to perform our simulations of the isolated H2O molecule on the MgO (001) surface. 

2.1 Calculation of the chemical potential using PMF 

The adsorbate molecules are all identical to each other, and are called species A; the atoms in the solid are called species B. 
We consider first the situation where A-molecules in the gas phase are in complete thermal equilibrium with those adsorbed on 
the surface. We assume that A-molecules cannot penetrate into the bulk. The situation we envisage is adsorption of molecules 
on the surface of a perfect crystal, with this surface corresponding to a particular crystal plane and being free of defects. This 
being so, we can draw a "reference plane" parallel to the surface, whose position is such that if translated outward by a few 
lattice spacings along the surface normal it would be entirely in the gas phase, and if translated inward by a similar amount it 
would be entirely in the bulk. The origin of coordinates lies in this reference plane, with the z-axis pointing along the outward 
normal. 
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When the entire system is in thermal equilibrium at a given temperature T and a given chemical potential of the A-molecules, 
there is a density po (number per unit volume) of A-molecules in the gas phase, and a surface density a (number per unit area) 
on the surface. In order to make a well defined, we have to decide what is meant by an A-molecule being "on the surface" . In 
order to do this, we define a "monitor" point for each A-molecule. If the molecule is an atom, this monitor point is just the 
position of the nucleus. If it is a molecule, it could be chosen as the position of a specified nucleus, or as the the centre of mass, 
or in some other way. The particular choice of monitor point will not (significantly) affect any of the physical results of interest. 
To separate A-molecules in the gas from those on the surface, we define a "transition plane", which is parallel to the reference 
plane, but displaced from it by the ^-displacement Zo- An A-molecule is counted as on the surface if the z-coordinate of its 
monitor point is less than zo, and in the gas phase otherwise. To ensure that the value of a does not depend significantly on 
the position of the transition plane, we choose zo so that the interaction of an A-molecule with the solid is negligible for z > zo. 

Since we have full thermodynamic equilibrium, the chemical potentials p sas (po,T) and p a ds(0", T) in the gas phase and on 
the surface must be equal. We are primarily interested in situations where po is so low that interactions between A-molecules 
in the gas are negligible. It is then convenient to express p ga s as: 

Mg as(po,T) = fc B Tln(poA 3 ) +mLs( t ) . ( 2 ) 

where the first term on the right is the chemical potential of a perfect gas of structureless particles, whose mass Ma is equal 
to that of the A-molecules, and A is the thermal wavelength A = {tf I'I-kMpJzb.T) 1 ^' 2 ■ If the A-molecules are atoms, then 
A i gas(7 1 ) = 0, but if they are molecules, Pgas^) represents the T-dependent contribution to p gaB (po,T) from internal vibrations 
and rotations. In a similar way, it is convenient to express p a ds(&, T) as: 

Mads (a, T) = k B T ln(aA 3 /d) + p[ ds (a, T) . (3) 

where d is an arbitrary fixed length. By analogy with the gas-phase chemical potential, the first term on the right represents the 
chemical potential of a gas of structureless particles of mass Ma confined to the surface region. The "excess term" p ads (c, T) 
depends, of course on the value chosen for d, but this is a convenient way to write Mads, because in the limit where quantum 
nuclear effects can be ignored, p^j^cr, T) is then independent of Planck's constant. This way of separating p ads is particularly 
helpful in the limit of low coverage, a — > 0, when /z ads becomes independent of a and includes the adsorption energy of A- 
molecules, as well as the entropy effects due to translations, vibrations and (hindered) rotations. For non-zero <r, /xL s also 
includes the energetic and entropic effects of adsorbate-adsorbate interactions. With these definitions, the combination of 
eqns @ and with the equilibrium condition p gas = Mads gives: 

cr/p = dexp [j3Ap 1 {a, T)] , (4) 

where Ap^ = p gaB — p ada - We refer to Ap^ as the "excess chemical potential difference" (ECPD). The adsorption isotherm (a 
as a function of po or gas pressure at constant T) can then be found by solving eqn Q for a. In calculating the thermodynamic 
properties of the adsorbate on the surface, it is therefore convenient to focus on the quantity ApJ (a, T) . 

We shall need later the relation between Apv and the isosteric heat of adsorption hi so , which can be defined as the negative 
slope of ln(po) plotted against 1/T at constant surface density a. From eqn Q, we have: 

^° = -(^HwW) =(^(/?V)) • (5) 

To develop a strategy for calculating Apy , we consider the spatial distribution of A-molecules in thermal equilibrium. Let 
p(r) dr represent the probability of finding the monitor point of any A-molecule in volume element dr at position r. If r n is the 
dynamical variable representing the position of the monitor point of A-molecule n and there are u A-molecules in the system, 
then p(r) is: 




., v ■ 

where (•) denotes thermal average. When r = (x,y,z) is near the surface, p(x,y,z) depends on x and y as well as z, and 
exhibits the translational periodicity of the surface. We are not interested here in the dependence on x and y, so we consider 
the distribution p(z), which is the spatial average of p(x,y,z) over x and y. For z far from the surface, p(z) is equal to the 
gas number density po, but in the region of the surface, p(z) will have a large peak (or perhaps more than one peak), due to 
A-molecules adsorbed on the surface. With the position zo of the "transition plane" chosen as above, this peak is entirely in the 
region z < zo, and for z > zo, p(z) is very close to po- The surface density a of adsorbed molecules is the integral of p(z) over 
the region z < zo. It is convenient to work with the distribution y(z) = p(z)/po, which is normalised so that lim z _,oo y(z) = 1. 
Then we have: 

a/p = / dzy(z) . (7) 



From this, we see that y(z) is closely related to the ECPD Apy . In fact, from eqn Q: 

exp [/3A M f (a,r)] = i f ° dzy(z) . (8) 
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Standard m.d. simulation can be used to calculate the unnormalised distribution p(z) = poy(z) in the region of z where the 
adsorbed molecules spend most of their time. Accumulation of a histogram for the probability distribution of z suffices for this 
purpose. However, in order to perform the integral in eqn (0 we need the normalised distribution y(z), and for this, the entire 
region z < zq must be sampled. Under most circumstances, simple accumulation of a histogram is not a practicable way of 
doing this, because the probability of adsorbed molecules sampling the region z ~ zo is so small. There are several well known 
techniques for overcoming this "rare-event" problem. The technique employed in the present work uses the potential of mean 
force (PMF) [17]. In this approach, y(z) is expressed in terms of a PMF <f>(z), which plays the role of a z-dependent free energy: 

y(z) = exp l-p<t>{z)] . (9) 

By standard arguments, the z-derivative d<j>/dz is equal to minus the thermal average of the z-component of the force acting 
on the monitor point of a chosen A- molecule. Denoting by (T z ) z the thermal average of the z-component of the force acting 
on the molecule when the z-component of its monitor point is constrained to have the value z, we then have: 

<P(z) = / dZ (T z ) z , , (10) 



where {T z ) z is counted positive if directed outwards. Combining eqns (JSJ and JSJ, we obtain a formula for the ECPD Ap/(cr, T) 
in terms of the mean force: 



A//(ff,T) = fc B Tln|iy ° dz e- Wz) | = k B Tlnl^J ° dz exp -0 dz'(T z ) z , j 



(11) 



The physical content of this expression is that the ECPD is determined by the reversible work (integral of mean force) performed 
on transporting the molecule from one phase to the other. The calculation strategy is to evaluate the mean force {T z ) z as 
a time average in a series of constrained DFT m.d. simulations at a sequence of z values, and to compute the mean force 
integral numerically. This strategy can in principle be applied at any temperature and coverage, provided ways can be found 
of achieving adequate statistical accuracy. It is one of the purposes of this paper to test the practical feasibility of the strategy. 

We note in passing an alternative way of calculating y(z), which should also be effective, namely umbrella sampling |18U19] . 
In this, the probability distribution of the dynamical variable being sampled is deliberately biased by adding to the Hamiltonian 
a potential that acts on this variable. In the present case, we would add a potential V np (zi) acting on the z-coordinate of the 
monitor point of a chosen molecule. It is an exact result of classical statistical mechanics that the probability distribution p(z) 
of the chosen molecule is related to the distribution y(z) in the absence of V up by p[z) = Ay(z) exp(— /3V up (z)), where A is a 
constant. If V up (z) is appropriately chosen (in particular, if it is similar to —<f>(z)), p(z) can be fully sampled over the required 
region z < zq by accumulation of a histogram, and y(z) can then be recovered by multiplying by exp(/3V up (z)). 



2.2 Thermal desorption rate 

The rate at which molecules desorb from the surface when there is complete thermal equilibrium between gas and adsorbed 
molecules can be derived by standard detailed-balance arguments, which we recall briefly. (The arguments are well known [20] . 
but we need to express them in a way that is consistent with the present notation.) According to elementary statistical 
mechanics, the outward flux k of molecules across the transition plane is k = po (k B T /2n Ma) 1 ^ 2 , where we make a negligible 
error in setting the number density at the transition plane equal to po- If the sticking coefficient S of molecules arriving from 
the gas phase is unity, then the flux n is entirely due to spontaneously desorbing molecules. But if S < 1, then only the flux 
Sk is due to spontaneous desorption. We denote by 7 the probability per unit time that a molecule spontaneously desorbs, so 
that 7<t = Sk. It follows from eqns @ and Q that: 

7 = {p /a)S(k B T/2nM A ) 1/2 = S (k B T/2nM A ) 1/2 / j ° dzy(z) 

I J — OO 

= S(k B T/27YM A ) 1/2 iexp [-/3A//] . (12) 

Although eqn (|12l) is derived for conditions of complete thermal equilibrium, it may sometimes be correct for the desorption 
rate in a TPD experiment, even though this is an irreversible process. A necessary condition for it to be correct is that the rate 
of equilibration of the adsorbate on the surface be fast compared with the desorption rate 7. In general, the desorption of a 
molecule from the surface leaves behind a distribution of the remaining molecules that is not typical of thermal equilibrium. If 
the molecules in the region where the desorption happened do not have time to equilibrate before the next desorption occurs 
in this region, then the use of eqn (|12|l is no longer strictly correct. The equation can also fail for another reason. In a TPD 
experiment, the molecules are deposited on the surface at a very low temperature, and T is steadily increased. If the rate of 
temperature increase is so fast that the adsorbate system has no time to equilibrate, then the desorption rate will be history 
dependent, and eqn (I12f) will fail, even at low coverage. In the simulations presented in Sec. [3j we shall study the relaxation 
rates associated with diffusion and reorientation of the H2O molecule on MgO (001), in order to test the correctness of the 
formula for 7. 
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2.3 Practical calculation of PMF 



It is clear from the foregoing theory that the primary quantity to be calculated is the mean force (J- Z ) z , from which we obtain 
the PMF and chemical potential by integration. The calculation of (T z ) z by m.d. simulation of H2O on MgO (001) raises a 
number of practical issues. We recall first that we are free to choose the molecular "monitor position" in different ways. In the 
present work, we choose it to be the position of the water O atom. At each time step of m.d., the electronic-structure calculation 
yields a Hellmann-Feynman force on the core of the water O exerted by the valence electrons and the other ionic cores. In 
constrained m.d. in which the ^-coordinate of water O is held fixed, T z is simply the z-component of this Hellmann-Feynman 
force. However, we note a disadvantage of this choice of monitor point, which is that T z does not vanish even when the molecule 
is far from the surface, because the water O oscillates with the internal vibrational modes of the molecule. This means that 
time averaging is needed to calculate {J- Z ) z , even when the constrained value of z is far from the surface. This problem can 
be avoided, and the statistical efficiency somewhat improved, by choosing the monitor point to be the centre of mass of the 
molecule, as we shall discuss in detail in paper II. 

Given our stated aim of calculating Afi^ correct to better than 50 meV for a chosen exchange-correlation functional, we 
need to give thought to the tolerance on the statistical error in the calculation of (J- Z ) z at each z- value, the number of z-values 
at which (T z ) z is needed, and the way in which the integrations of eqns (|10|1 and should be performed. Tests show that in 
the absence of statistical errors the integration of eqn (|10|l can be performed to obtain 4>{z) over the relevant z-range to much 
better than the required tolerance if we have {T z ) z at ~ 15 roughly equally spaced values, and if the integration is performed 
by the trapezoidal rule. In the integral f*° dz exp(— f3<j){z)) needed to obtain A//, the integrand varies rapidly in the region of 
its maximum, but we obtain the required accuracy by performing a cubic spline fit to the the <f>(z) values. Coming to statistical 
errors, we derive in Appendix A a relationship between the errors on the (T z ) z values and the resulting statistical error on A//, 
which allows us to estimate in advance the required duration of the constrained m.d. simulations needed to reduce the error 
below our required tolerance. 

The m.d. simulations were performed in the canonical (N,V,T) ensemble, with a Nose thermostat [3T]. However, careful 
attention needs to be give to the issue of ergodicity, because the MgO substrate is close to being harmonic at all temperatures of 
interest here, and the internal H2O vibrations are not only nearly harmonic but also of much higher frequency than the lattice 
modes. To ensure efficient transfer of energy between all degrees of freedom, we add to the Nose thermostat also an Andersen 
thermostat 22 , in which the velocities are periodically randomised by drawing new velocities from the Maxwellian distribution 
appropriate to the desired temperature. 

In considering how to achieve the required statistical accuracy on {T z ) z at each z, it is important to note that, since {T z ) z 
is a static thermal average, it is completely independent of the atomic masses. This means that the masses used to generate 
the m.d. trajectories do not have to be set equal to their physical values, but can be chosen to improve the efficiency with 
which phase space is explored. This is important for the present system, because of the very high vibrational frequencies of the 
H2O molecule. By artificially increasing the H mass, we can take a larger m.d. time step, without affecting the final results. A 
convenient way to gauge the advantage gained by altering the H masses is to consider the number of m.d. time steps needed 
to achieve a specified statistical accuracy in the calculation of (T z ) z for a given value of z, the statistical accuracy being given 
by the re-blocking analysis described in Appendix A. In test calculations, we found that if mH is increased from 1 to 4, the 
time step can be increased from 0.5 to 1.0 fs, and the number of time steps needed to achieve the same accuracy is reduced by 
a factor of two. Further increase of win makes little difference to the sampling efficiency. The m.d. simulations reported here 
were performed with the choice mu = 8, the masses of O and Mg atoms being their physical values. 

Our ab initio m.d. simulations were performed with the projector-augmented- wave (PAW) implementation of DFT |23l I24j . 
using the VASP code [55]. The plane- wave cut-off was 400 eV, and the augmentation-charge cut-off was 605 eV. We used core 
radii of 1.06 A for Mg and 0.80 A for O. There have been many previous DFT studies of H 2 on MgO (001) [26], and it is 
generally agreed that the molecule does not dissociate at low coverage. The interaction of the molecule with the surface is partly 
electrostatic, but there is an important contribution from non-bonded interaction of the O lone pair with surface Mg ions, and 
hydrogen bonding of H with surface O. The choice of exchange-correlation functional for such non-bonded interactions is not 
a trivial matter, and it can make a large difference to interaction energies, as will be shown in the next Section. Most of our 
calculations are performed with the Perdew-Burke-Ernzerhof (PBE) functional [27], which is usually taken to be one of the best 
parameter-free forms of exchange-correlation functional. The simulations employ the usual slab geometry, and the simulation 
conditions were decided on the basis of preliminary tests, which are described next. 

3 Simulations: H2O on MgO (001) at low coverage 
3.1 Preliminary tests 

Since our aim is to calculate the chemical potential and desorption rate of an H2O molecule in the limit of zero coverage on the 
surface of semi-infinite bulk MgO, we have made tests to ensure that the slab thickness, vacuum gap and size of the surface 
unit cell are all large enough to bring the system very close to this limit. The tests were done on the adsorption energy -E a ds, 
defined in the usual way as the sum of the energy _E(H20) of an isolated water molecule and the energy _E(MgO) of the relaxed 
clean MgO slab minus the energy _B(H20 + MgO) of the relaxed system in which the molecule is adsorbed on the surface of the 
slab: 

£ ads = E(K 2 0) + £(MgO) - E(R 2 + MgO) . (13) 
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With this definition, _E a ds is positive if the total energy decreases when the molecule goes from the gas phase to the adsorbed 
state. 

With PBE calculations, the most stable adsorbed configuration we have found (Fig.QJ,,b) has the water O bound to a surface 
Mg, one of the two O-H bonds directed to a surface O, and the other tilted up at an angle to the surface plane. In this 'tilted' 
geometry, -E a< j s is converged to within ~ 1 meV with a slab containing three layers of ions and vacuum width of 12.7 A (this is 
the distance between the layers of surface ions that face each other across the vacuum gap). With a 2 x 2 surface unit cell (16 
ions per layer in the repeating cell), -E a ds is converged to better than 5 meV. With these settings, F-point sampling achieves 
convergence to within 5 meV. The static adsorption energy given by PBE in this geometry is 0.46 eV. The calculations also 
reveal a second stable adsorbed configuration having a slightly lower E^ of 0.45 eV, in which the water O is bound, as before, 
to a surface Mg ion, but with the molecular plane almost flat on the surface, and with both O-H bonds directed to surface 
O ions (Fig. QJ,d). When the calculations are repeated with the LDA, we find that for the 'flat' adsorbed geometry -E a d s is 
0.95 eV. The large difference between the PBE and LDA values of J5 a d s indicates that quantitative agreement between theory 
and experiment cannot be expected without calibrating DFT calculations against more reliable methods. 

Guided by these tests, we have performed all the following calculations with the 3-layer 2x2 slab (48 ions in the MgO slab 
per repeating unit), a vacuum gap of 12.7 A, and T-point sampling. Only the PBE approximation is used from now on. At all 
temperatures, we set the lattice parameter equal to 4.23 A, which is the T = K value in the bulk crystal, according to PBE. 

3.2 Results for PMF, chemical potential and desorption rate 

We have stressed that the feasibility of the present scheme depends on being able to make the statistical errors on (T z ) z small 
enough with simulation runs of affordable length. We use the analysis presented in Appendix A to determine the statistical 
errors on A/r. We find that at 400 K, with typically 15 z-points and with runs of 40 ps duration at each z-point, the statistical 
error on A^ is less than 20 meV, which is much better than the accuracy we are aiming for. The statistical errors are similar 
at other temperatures. 

The mean force and the PMF from our simulations at four temperatures are reported in Figs. [2] and [3] We note the strong 
temperature dependence of both quantities. Interestingly, (T z ) z shows a double minimum at low T, which we believe is due 
to the fact that the characteristic orientation of the H2O molecule changes substantially with z. The temperature dependence 
of the well-depth in <\>{z) is closely related to the entropy part of A/J , as can be seen by representing <f>(z) in the simple form 
4>(z) = m in + |c(« — Zmin) 2 , where m in is the (negative) value of <j>(z) and c is the curvature at the bottom of the well. With 
this approximation, we have: 

A/i f = fc B Tln [d- 1 (27rfc B r/c) 1/2 exp(-/30 min )] = -0 min (T) + h B Tln(2nk B T/cd 2 ) . (14) 

Then the entropy s — —dA/i'/dT is: 

s = --ka (l + \n(2nk B T/cd 2 )) + d<f> mla /dT . (15) 

The first term on the right represents the entropy contribution from confinement of the z-coordinate of the molecule. The 
second term, which is positive (the well-depth decreases with increasing T), is due to the confinement of its translational and 
rotational degrees of freedom, as will be discussed in more detail below. Our calculated values of A/j} at the four simulation 
temperatures vary rather linearly with T, and this implies that the isosteric heat of adsorption hi BO deviates only slightly from 
the T = K value of the adsorption energy -E a d s - 

In order to obtain the temperature-dependent desorption rate 7, we need to know the sticking coefficient S. We expect this 
to be close to unity, because the PMF shows that there is no barrier to adsorption. To test whether S deviates significantly 
from zero, we have conducted a series of simulations in which the 3-layer MgO slab has already been well equilibrated at a 
given temperature T; the H2O molecule is initially in the middle of the vacuum gap, and its atoms are given random velocities 
corresponding to that T. The subsequent time evolution is then monitored. We show in Fig. Utile results of 12 such simulations 
at T = 400 K. We see that in all cases H2O becomes adsorbed on the surface, and there is no indication of desorption in the 
period of several ps after adsorption, during which time the molecule will have become well equilibrated on the surface. In the 
light of this, it is an accurate approximation to set 5=1. 

Numerical values of the desorption rate 7 as a function of T calculated from eqn (I12|l are shown as an Arrhenius plot in 
Fig. [5] We have only a few data points, but the Arrhenius plot appears to rather straight, except at the highest T. This 
means that our results are fairly well described by the Polanyi-Wigner formula, eqn fTJ. We note that for the present case of 
low-coverage desorption in the absence of dissociation, the reaction order n is unity, so that the Polanyi-Wigner formula becomes 
7 = —6~ 1 d9/dt = / exp(— AE/k B T). If we pass a straight line through the points at the two lowest temperatures in fig. [5] we 
obtain an activation energy AE = 0.454 eV, and a prefactor / = 2.7 x 10 15 s _1 . We note that the activation energy is almost 
exactly equal to the zero-temperature static adsorption energy £? a ds = 0.46 eV. In the temperature range 200 — 300 K, where 
thermal desorption of H2O from MgO (001) is experimentally observed at sub-monolayer coverage, the calculated prefactor of 
2.7 x 10 15 s _1 is thus enhanced by a factor of ~ 270 above the value of 10 13 s _1 often used in analysing TPD measurements on 
this system in the past |281 129] . 

The enhancement of the prefactor arises from the T-dependence of the PMF. On substituting the approximate form of AfJ 
given by eqn (|14|) into eqn (|12p . we find: 

7 = S (c /4tt 2 M a ) 1/2 exp(/30 mi n) . (16) 
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The factor (c/47t 2 Ma) 1 ^ 2 is the oscillation frequency of the centre of mass of the molecule along the surface normal. If (pmiu varies 
roughly linearly with T, then ^ min = -_E ads + T(d<p min /dT), so that 7 = S{c/An 2 M A ) 1/2 exp(fc- 1 ^ m in/9T) exp(-/3£ ads ). The 
enhancement factor is thus exp(fcg 1 (90 m i n /9T), and it comes from the part of the adsorption entropy associated with confinement 
of the translational and rotational degrees of freedom other than oscillation along the surface normal. 

3.3 Spatial distribution and memory 

We now discuss two important questions. The first concerns the more detailed interpretation of the prefactor enhancement. 
Since this is due to confinement of degrees of freedom of the molecule, it should be possible to relate it to the probability 
distributions of these degrees of freedom. The second question concerns the crucial assumption underlying our calculation of 
the desorption rate, that equilibration rates of the adsorbed molecule are faster than the desorption rate. Equilibration rates 
are related to memory times, and we want to use our simulations to characterise these memory times. 

Since H2O contains three atoms, it has nine degrees of freedom. One of these is bodily vibration along the surface normal, 
whose distribution is y{z). The three internal vibrational modes of the molecule do not undergo large changes on adsorption 
at low coverage, so that they do not contribute significantly to the enhancement of /. The enhancement therefore comes from 
the remaining degrees of freedom, two of which refer to translation in the surface plane, and the other three to rotations. We 
describe in Appendix B how the enhancement factor can be related semi-quantitively to the non-uniformity of the probability 
distributions of these degrees of freedom. 

We examine the reduced freedom of the translational motion by plotting the probability distribution pt T (x,y) of the x- 
and ^-coordinates of the water O atom. This is done by dividing the surface plane into a square grid, and accumulating the 
frequencies with which x and y fall in each cell of the grid. Since p tT (x,y) has the translational symmetry of the crystal surface, 
we can improve the statistics by averaging over symmetry related cells. From our calculated pti{x,y) at T = 400 (Fig. EJ), 
we see that the water O spends most of its time close to surface Mg sites, in line with our finding (Sec. 13. ip that the most 
stable relaxed configuration of H2O on the surface has water O coordinated to surface Mg. However, it is also clear that the 
most probable O position is not directly above Mg, but is displaced to positions roughly along a cubic axis in the surface 
plane. This agrees with our static calculations on the tilted configuration.. The analysis of Appendix B shows that we can 
estimate the entropy reduction associated with p tI (x,y) as follows. We fit the calculated pti(x,y) with a superposition of four 
Gaussians A exp (— a|r — R| 2 ) , centred at the positions Sx(±l, 0) and 5x(0, ±1) relative to the Mg site, using a and Sx as fitting 
parameters. We then normalise the resulting smoothed distribution so that J cell dx dy p tv (x,y) — 1 where the integral goes over 
the surface unit cell, whose area is A cc n. Denoting by ci the constant such that the maximum value of cipn(x,y) in the cell 
is unity, the contribution to the enhancement factor from translational confinement is estimated as bi = A ce \\/cy. At T = 300 
and 400 K, we find 61 = 3.8 and 2.9, respectively. 

The orientation of the molecule is specified by three angles 9, <p and ip. These characterise the orientation of the molecule 
relative to a reference orientation, which we choose to be such that the bisector of the H-O-H bond points along the outward 
normal, with the plane of the molecule lying in the x-z plane (the x-axis is along one of the cubic axes of the substrate). From 
this starting point, we can produce any other orientation in three steps: (i) rotate the molecule about the bisector by angle <f>; 
(ii) rotate the molecule by angle 9 about the axis perpendicular to the bisector in the molecular plane (the result of this is that 
the bisector makes an angle 9 with the surface normal); (iii) rotate the molecule about the bisector by angle ip. The angles 
(9, <j>, ip) are thus the conventional Euler angles specifying rigid-body rotation of the molecule. 

The probability distribution pbis(#, <p) of the bisector is computed by dividing the 9 range (0, n) and the <p range (0, 2n) into 
a uniform grid and accumulating a two-dimensional histogram. The contour plot of Pbis(9, 0) at T = 400 K (Fig. [7]) shows that 
the direction of the bisector is always quite close to the surface plane, and within this plane tends to point along one of the 
diagonal directions. (This is exactly the orientation of the 'flat' configuration (Fig.[T}, in which the two O-H bonds are directed 
from a surface Mg site to two of the nearest-neighbour O sites.) This being so, the angle ip effectively measures the tilt of the 
molecular plane relative to the surface. The rather broad calculated probability distribution p sp (ip) of the angle ip (Fig. [8}, 
centred on ip = 0, means that the molecular plane tilts relative to the surface over a rather wide range. The double-peaked 
structure of the distribution indicates that the untilted configuration is not the most stable, again in accord with our static 
results. 

The reduction of rotational freedom associated with (9, <p) is estimated by fitting the distribution pbis(#, </>) to a superposition 
of Gaussians. As explained in Appendix B, the contribution &2 to the enhancement factor due to confinement of 9 and (p is 
estimated by a procedure similar to the one used for translational confinement. At T = 300 and 400 K, we find 62 = 19.9 and 
16.0, respectively. Finally, analysing the confinement of ip shown in Fig. [8j we obtain enhancement contribution 63 = 2.8 and 

2.4 at T = 300 and 400 K, respectively. Combining all the translational and rotational factors, we arrive at total enhancement 
factor of 212 and 111 at T = 300 and 400 K, which are roughly consistent with what we obtained from our numerical results 
for 7. 

Turning now to the question of memory times, we want to give evidence about the dynamics of the two processes that 
are likely to have the longest memory times, namely hopping between different surface sites, and transitions between different 
orientations. We have attempted to calculate appropriate correlation functions, but with only a single molecule in the system, 
it is not possible yet to achieve good statistics. The data we present are therefore only semi-quantitative. 

To illustrate the translational dynamics, we show in Fig. [9] the time-dependent x(t) and y(t) coordinates of the water O 
over the 100 ps span of an unconstrained simulation at 400 K. We see several well-defined intersite jumps in which x and/or y 
change by the Mg-O nearest-neighbour distance d = 2.12 A. In most of these jumps, both x and y change, and this indicates 
a diagonal jump between Mg sites, but there are events in which one coordinate changes by 2d, with no change in the other. 
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By simple counting of jumps, we estimate a hopping rate of 1.4 x 10 11 s _1 . The same procedure at T = 300 K gives a hopping 
rate of 3.8 x 10 10 s _1 . To interprete these results further, we have made nudged elastic band calculations [3D] of the energy 
barrier for intersite hopping, and we find the value 0.13 eV. This is roughly consistent with the ratio of hopping rates of about 
3.7 between 300 and 400 K. It is clear from this that, in this temperature region, intersite hopping is very much more rapid 
than the desorption rate, and this will become even more true at lower temperatures. 

To display the rotational dynamics, we show in Fig. 1101 plots of the angles 6(t), <j>(t) and ip(t) during the course of the 
simulation at 400 K. The dynamics of 9 and if) consists of very rapid fluctuations on a time scale of 1 ps or less, over the 
range expected from the probability distributions of Figs. [7] and [5] The angle (j> has somewhat slower dynamics associated with 
hopping between the four equivalent sub-sites around each Mg site (Figs [1] and [TJl , but at T = 400 K, this hopping rate is 
~ 5 x 10 11 s~\ which is several times faster than the hopping rate between sites. 

The conclusion from this is that all equilibration rates on the surface appear to be very much faster than the desorption 
rate at temperatures of interest. 

4 Discussion 

The practical calculations we have presented show the feasibility of calculating thermodynamic properties of a surface adsorbatc 
by ab initio methods, without recourse to a lattice-gas approximation. We have assumed nothing at all about the adsorption 
sites of H2O on MgO (001), and the ab initio m.d. simulations themselves automatically sample the sites and orientiations 
that are statistically significant. In the same spirit, we completely avoid approximations such as the harmonic approximation, 
and we fully include the coupling of the molecular and substrate degrees of freedom. In this sense, there are no statistical- 
mechanical approximations whatever, except for the neglect of quantum nuclear effects, to which we return below. The single 
uncontrollable approximation is the DFT exchange-correlation functional. Furthermore, when calculating the desorption rate, 
we do not assume the validity of the Polanyi-Wigner formula, but we use this formula only as a means of fitting the computed 
results. In this way, we obtain an ab initio value of the frequency prefactor. 

It is well known from both experiments |31| and simulations |19] on a wide range of systems that frequency prefactors 
for thermal desorption often differ, sometimes by many orders of magnitude, from the value of 10 13 s _1 that might naively be 
expected if the prefactor is thought of as an "attempt frequency" . It is also well known that this is due to the strong reduction of 
translational, rotational and conformational entropy that often occurs when a molecule goes from the gas phase to the adsorbed 
state. For H2O on MgO (001), our calculated prefactor has the value / = 2.7 x 10 15 s -1 , so that it is enhanced by a factor 
of over 100 above the typical vibrational frequency of the molecule relative to the surface. We have seen that an enhancement 
factor of this general size is expected from the translational and rotational probability distribution functions of the molecule. 

A crucial assumption behind the Polanyi-Wigner formula is that the desorption rate depends only on the instantaneous 
temperature and coverage, and is not history dependent. A key condition for this to be true is that the equilibration rate of the 
molecule on the surface be faster than the desorption rate. We have attempted to characterise the typical equilibration times 
for the isolated H2O molecule on MgO (001) by studying the diffusional and rotational dynamics, and we have found that the 
condition is satisfied by a margin of several orders of magnitude in the temperature region where thermal desorption can be 
experimentally observed. 

A direct comparison with experimental TPD data for H2O on MgO (001) is not straightforward. There have been several 
TPD studies reported 28, 29 , some of which refer to carefully prepared surfaces that appear to be relatively free of defects. 
Desorption at sub-monolayer coverage is associated with a TPD peak at ~ 245 K. In our simulated system, the desorption rate 
at this T is ~ 2.6 x 10 6 s — l , so that all molecules would desorb in less than 10 jits, which is many orders of magnitude less than 
the time-scale of a TPD experiment. This might suggest that the adsorption energy of 0.46 eV given by the PBE functional we 
have used is considerably too low. However, the effect of attractive water-water interactions may be important even at coverages 
well below the monolayer level. In fact, adsorption isotherm measurements indicate a critical point in the surface phase diagram 
of H2O on MgO (001) at T ~ 210 K |32| . The effect of water-water interactions on the desorption rate clearly needs to be 
quantified. A further complication is that we have so far ignored quantum nuclear effects. Because of the very high vibrational 
frequencies of the water molecule, it is possible that changes of zero-point energy on adsorption might shift i5 ac i s significantly. 
Even without all these effects, it is not clear that we should expect good agreement with experiment yet, because the calculated 
Ends depends so much on the exchange-correlation functional. The frequency prefactor / might also depend significantly on 
exchange-correlation functional. If we assume provisionally that our calculated / of 2.7 x 10 15 s _1 is essentially correct, and 
we ignore water-water interactions, then the experimental TPD peak temperature of 245 K would require an activation energy 
AE = 0.78 eV, which is well above the PBE adsorption energy of 0.46 eV, though it is still below the LDA value of 0.95 eV. It 
is clear that DFT predictions of I? a ds need to be tested against more accurate and reliable methods. 

A number of future challenges are suggested by this work. The most obvious of these is the extension of the calculations 
to higher coverages. According to the theory we have presented, the calculation of the PMF on a chosen molecule, and the 
integration of the resulting distribution function y(z), allows us to calculate the chemical potential and the desorption rate at 
arbitrary temperature and coverage. At the time of writing, we have performed exploratory ab initio calculations of this kind 
for H2O on MgO (001) at coverages of 0.25 and 0.5 ML. However, the problem that emerges is that the memory times are 
much longer than for the isolated molecule, so that considerably longer simulations are needed in order to achieve acceptable 
statistical accuracy. This problem of statistical sampling becomes rapidly worse at low temperatures, and so far we have 
achieved stastically accurate results only at T = 800 K, which is well above the region of practical interest. A second important 
challenge is that of going beyond DFT. We have noted that the LDA and GGA forms of exchange-correlation functional give 
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static adsorption energies differing by roughly a factor of two. This means that we must envisage future calculations of the 
present kind, but performed by post-DFT techniques. The recent proposal of a method for performing quantum Monte Carlo 
simulations along an m.d. trajectory |33| may indicate one way to do this. Recent successes in applying high-level quantum 
chemistry to condensed-matter energetics are also promising [34]. A third challenge is that of eliminating the approximation 
of classical statistical mechanics for the nuclei. For the case of H2O on MgO (001), the errors due to the use of the classical 
approximation will not be large, but may be significant if one wishes to achieve chemical accuracy. We will report in the 
second paper of this series on the generalisation of the present theory to quantum statistical mechanics for the nuclei, using 
path-integral ab initio simulation. 



5 Conclusions 

In summary, we have shown how ab initio methods can be used to calculate the chemical potential of an adsorbate, with 
full inclusion of entropy effects. The methods used can in principle be applied at any coverage. We have shown how the 
methods also yield values for the desorption rate and hence the frequency pre-factor in the Polanyi-Wigner formula. For the 
case of H2O on MgO (001) at low coverage, this pre-factor is enhanced by at least two orders of magnitude above the values 
generally assumed in the past, and we have given a detailed interpretation of this enhancement in terms of the confinement of 
translational and rotational degrees of freedom. The crucial condition of rapid equilibration necessary for the validity of the 
Polanyi-Wigner formula appears to be satisfied by a wide margin for low-coverage H2O on MgO (001). Preliminary comparisons 
with experimental data suggest that for this system the adsorption energy given by PBE may significantly too low. 



Appendix A: Analysis of statistical errors 

We explained in Sees. I2.3l and l3.2l the importance of monitoring the statistical errors in the calculation of the mean force {J- Z ) z 
and hence of the PMF and the chemical potential. We summarise here how we have estimated the errors on the mean force, 
and how we have used this to estimate the errors on the other quantities. 

The mean force is calculated at a set of z- values: z.o > zi, . . . > z n - At a given z- value Zi, the estimate Ti of the mean 
force obtained by averaging over the length of the run differs from the exact value Tf^ that would be obtained if the sampling 
were perfect. We denote by 8Ti = T% — Tf^ the difference that occurs in a given simulation run. We estimate the standard 
deviation (STf) 1 ^ 2 by the usual re-blocking method. In this method the simulation of total duration r is divided into v 
blocks, each of duration t/v, and we compute an estimate of Ti for each block I. Then, we compute the quantity 

a 2 . — v~ x 2^=1 0^(0 ~~ ^) ' where Pi i s the estimate obtained by averaging over the entire time r. If the duration t/v of 
each block is short, the block averages axe strongly correlated, and <t„ underestimates the true statistical error. However, 
as the duration of the blocks becomes longer than the correlation time, the become statistically independent, and o v tends 
to a plateau value. The reblocking technique consists of plotting o v against t/v and taking the standard deviation (ST 2 ) 1 ^ 2 to 
be the plateau value of o v as v is decreased. 

We now turn to the statistical error on the PMF, denoting the value of <j>(z) at the ith z-point obtained from a given set of 
simulation runs by <j>i, the exact value by and the difference <j>i — (jyf 1 by 8<f>i. Since the <f>i values are obtained by integrating 
inward from the largest z-value zq using the trapezoidal rule, we have: 



2 

n=l 

1 1 * _1 1 

bo + -(zo - z\)To + g / - z j+1 )Tj + ^{zi-i - zi)Ti . (17) 



3=1 

The errors on Tj and jF fc are statistically independent for j 7^ k, so that: 



-1 



(zo — zi) (SFq) + y y (zj-i — Zj+i) (8J r j ) + (z i - 1 -z i ) (5Ti) 
2=1 



(18) 



Now to estimate the statistical error on , we recall (eqn (|11[) 1 that = kBTlnY, where Y = d 1 dz exp(— j34>{z)). 
In practice, the integral comes almost entirely from a narrow region of z around the value Zmin where 4>(z) has its minimum. But 
the statistical fluctuations of <j>i at Zi that are near each are almost perfectly correlated, since the fluctuation 5<j>i comes from 
the accumulation of fluctuations 8Tj at many points Zj > Zi. So to calculate the fluctuation of Y, it is a good approximation to 
say that the fluctuations of all <f>i in the region of z mln are perfectly correlated. This means that 5Y — — P5<f)minY , where <5<^ m m 
is the fluctuation of S<j>i for Zi — z m in- Then since SAfi^ = ksTSY/Y, we finally obtain the estimate for the standard deviation 
((5A^) 2 ) 1/2 = (<5<#L ln ) 1/2 , with ((5^ in ) 1/2 calculated according to eqn {TSJ). 
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Appendix B: Trans lational and rotational distributions and enhancement 
of the frequency prefactor 

We noted in the text that the entropy of the adsorbed molecule is less than that of the free molecule, because its translational 
and rotational freedom is reduced. This entropy reduction is closely related to the enhancement of the frequency prefactor / 
above its naively expected value of ~ 10 13 s _1 . We want to make a semi-quantitative connection between the entropy reduction 
and the translational and rotational distributions presented in Sec. 13.31 

For this interpretative purpose, we ignore the degrees of freedom of the substrate, and we assume that the x-y translational 
distribution, the ^-distribution, the 9-<f> distribution and the ip distribution are all independent of each other. This is equivalent 
to supposing that the distribution over x, y, z, 9, (j> and tp is governed by a Boltzmann factor exp(— (3V(x, y, z, 9, <f>, ip)), where 
V is expressed as: 

V(x, y, z, 9, </>, VO = -Vo + ^az 2 + u(x, y) + v{9, 4>) + w (■>/)) . (19) 

Here, — Vo represents the energy of the molecule in its most stable adsorbed configuration, which implies that u > 0, v > and 
w > 0. Then the distribution y(z) defined in Sec. 12. H is given by: 



y(z) = Cexp 



-0 (-Vo + \az' 



(20) 



■>7T 



x / dxdye- 0u{x ' v) I d<j> I d9 sm9 e - 0v{eM / d^e~ l3wW (21) 

./cell JO JO Jo 

for z near the surface, where the x-y integration goes over a single surface unit cell. By definition, y(z) = 1 far from the surface, 
and this fixes the constant C: 

1 = 87r 2 CA ccll , (22) 

where A ce n is the area of the surface unit cell. In this approximation, the factor B by which the frequency prefactor is enhanced 
is the product of three factors: B = 6162^3 , where: 

61 = Accii / / dxdye-^'^ 

I J cell 

//■2-7T ("IX 
J d(f>jd9 sinfle-^ 8 '^ 

63 = 27V / J d ^ e ~ 0WW ■ ( 23 ) 
The potential of mean force <f>(z), defined by y(z) = exp(— f3cj>(z)), is then given by: 

4>{z) = -Vo + ^az 2 + k B T ln(6ife 2 6 3 ) , (24) 

so that the value of z at the minimum of <j) and the curvature at the minimum remain the same, but the well becomes less deep 
with increasing temperature, because of the factor fcBTln(&i&2&3). 

We estimate the values of the factors bi from the probability distributions as follows. We calculate the probability distribution 
Ptr(x, y) by histogram accumulation, as described in Sec. 13.31 and normalise it so that J" dx dyptr(x, y) = 1. We then determine 
the constant ci such that the maximum value of ciptr(x,y) in the cell is unity (ci thus has dimensions of area). We then have 
bi — A ce i\/c\. Similarly, with probability distribution pbis(#,<?i) normalised so that J d(j>d9 sin 9p^ s {9, <f>) — 1, we find the 
constant C2 such that the maximum value of C2Pbis(#, 4>) is unity. Then 62 = 47r/c2. Similarly for 63. 
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Figure 1: Top and side views of the 'tilted' ( panels a, b) and 'flat' (panels c, d) configurations of the H2O molecule adsorbed on 
the MgO (001) surface. Oxygen: large white spheres; Mg: medium grey spheres; H: small black spheres. 
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Figure 2: The mean force {T z ) z on the water O atom as function of its height z above the surface at T = 100, 300, 600 and 800 K 
(solid, dotted, chain and dashed curves, respectively, are guides to the eye). Bars on data points show statistical errors. Height 
z is relative to a fixed atom in the centre of the slab. 
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Figure 4: Set of 12 trajectories from simulations used to determine sticking coefficient S. Plots show ^-coordinate (A units) of 
water O atom relative to centre of vacuum gap between slabs. 
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Figure 5: Arrhcnius plot of the desorption rate 7 of H 2 from MgO (001) calculated using the PBE exchange-correlation 
functional. Bars on calculated values show statistical errors. The straight line is drawn to pass through the calculated values at 
the two lowest temperatures. 
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Figure 6: Contour plot of the spatial probability distribution of water O atom in the x-y plane at T = 400 K. Bottom right and 
top left corners of plot are Mg sites; bottom left and top right corners are O sites. Probability density is in arbitary units, with 
equal spacing between contours. 
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Figure 7: Contour plot of probability distribution of angles 9 and </> specifying orientiation of the bisector of H2O molecule (see 
text) at T — 400 K. Peaks of the distribution correspond to the four equivalent orientations in which the bisector is nearly parallel 
to the surface (9 ~ ^ir), and points along one of the diagonal directions (0 = jTt, jir, jir, and jir). 
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Figure 8: Contour plot of probability distribution of angle ip specifying rotation of H2O molecule about its bisector (see text) at 
T = 400 K. When the molecular bisector is parallel to the surface, the molecular plane is parallel to the surface when tjj — ^ir. 
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Figure 9: Time variation of the x- and y-coordinates of the water O atom in the course of an m.d. simulation at T = 400 K. 
Horizontal lines mark x- and y-coordinates of perfect-lattice sites, so that spacing between neighbouring lines is ^cjq = 2.115 A. 
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